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Abstract 

The present paper introduces a linear reformulation of the Kuramoto model describing a self- 
synchronizing phase transition in a system of globally coupled oscillators that in general have different 
characteristic frequencies. The reformulated model provides an alternative coherent framework through 
which one can analytically tackle synchronization problems that are not amenable to the original Kuramoto 
analysis. It allows one to solve explicitly for the synchronization order parameter and the critical point 
of 1) the full phase-locking transition for a system with a finite number of oscillators (unlike the original 
Kuramoto model, which is solvable implicitly only in the mean-field limit) and 2) a new class of continuum 
systems. It also makes it possible to probe the system's dynamics as it moves towards a steady state. While 
discussion in this paper is restricted to systems with global coupling, the new formalism introduced by the 
linear reformulation also lends itself to solving systems that exhibit local or asymmetric coupling. 
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I. INTRODUCTION 



Spontaneous synchronization of coupled oscillators with different natural frequencies is at the 
core of many striking phenomena in dynamical systems in realms from biology and physics to 
social dynamics The complexity of these systems defied attempts to encapsulate them in 

tractable mathematical formulations until, in 1975, Kuramoto produced a model that he was able 
to solve exactly for systems containing an infinite number of globally weakly coupled nonlinear 



oscillators Q3|, |4|, |5L |6fl. In such a system, the global coupling strength across the oscillators and 
the width of the oscillators' initial characteristic frequency distribution determine whether or not 
the system will self- synchronize, and with these Kuramoto was able to describe a phase transition 
into a self-synchronizing system. The Kuramoto model has since provided the basis for many later 
efforts to explore spontaneous synchronization. More importantly still, it is a valuable explanatory 
model in its own right for understanding numerous situations involving synchronization. In addi- 
tion to examples in biology (see IllE] and reference therein) such as neural firing patterns, clouds 
of fireflies flashing as one, and the coordinated action of cardiac pacemaker cells, many examples 
of the Kuramoto model have been recently discovered throughout the physical sciences including 
in such diverse systems as Josephson junction arrays |7], collective atomic recoil lasing [8], fla- 
vor evolution of oscillating neutrinos [9], and phase locking of oscillations in coupled chemical 
reactions lioll . 

This paper reformulates the Kuramoto model in terms of linear dynamics, permitting its so- 
lution through an eigenvalue/eigenvector approach. The analysis here is restricted to solving for 
the critical point of the fully locking transition and the synchronization order parameter beyond 
this transition; the properties of partially locked states are not considered. Within this regime, the 
reformulation of the Kuramoto model has a number of advantages over the original model. In 
particular, in the continuum limit, this linear model makes it possible to solve exactly a new class 
of synchronization models that is distinct from the class of systems for which the Kuramoto model 
in its original form has an analytic solution. Furthermore, whereas in the original version of the 
Kuramoto model the synchronization order parameter is only solvable in the continuum limit and 
then only implicitly, the alternative form presented below allows the order parameter to be solved 
explicitly and for any number of oscillators. In addition, the linearity of the reformulation makes 
it possible to investigate the time evolution of a system's self- synchronization and lends itself to 
adaptation to systems that exhibit local and/or asymmetric coupling between oscillators. 
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The present paper will begin with an introduction to the original Kuramoto model and its im- 
plicit analytic solution in the thermodynamic limit. Then the linear reformulation and its solution 
will be presented, and the mapping between the Kuramoto model in its currently accepted form 
and linear version presented here will be shown. This will be followed by some examples to 
demonstrate properties of the linear reformulation. 



H. KURAMOTO MODEL 

The Kuramoto model [3,4,5,3] describes a collection of N oscillators that are weakly coupled: 

JV 

flfc = u h + ^ Kjk sin(0j - 6 h ) (1) 

where Kj k is the coupling constant, uj k is the characteristic frequency of the k th oscillator, and 9 k 
is its phase. Although Kuramoto's method only yields an analytic solution for a uniform coupling 
scheme (further discussion below), note that the model's coupling constant is general and does 
not imply uniform coupling. Kuramoto was able to exactly solve eq. (OQ) for a system of globally 
and uniformly coupled oscillators, i.e. Kj k = K/N, with the constraints that the system be in the 
thermodynamic limit, i.e. N — > oo, and have reached a steady-state; in so doing, he showed the 
existence of a nontrivial self-synchronization phase transition in such a system. In this paper, as a 
measure of the synchronization of a system at a point in time, we will use the amplitude r of the 
complex phase order parameter, 



(where </> represents the collective phase of the synchronized state) which ranges between (no 
synchronization) and 1 (perfect synchronization); steady-state synchronization of the system will 
be specified as such. Kuramoto's result was an implicit equation for r in the synchronized state 
given by 

/ZL 
2 cos 2 6g(Kr sin 6)d6 (3) 
■f 

where g = g(u> k ) is the distribution of co k and the frame of reference is the rotating frame in which 
the frequency of the synchronized solution is zero. We will refer hereafter to eq. ©, with the 
abovementioned accompanying constraints, as the Kuramoto solution. The Kuramoto model eq. 
O, along with its many variations, has been studied in great detail |fj], and is seen as the standard 
model for synchronization of globally weakly coupled nonlinear oscillators. 



III. LINEAR REFORMULATION AND SOLUTION 



We now propose a linear reformulation of the Kuramoto model of spontaneous synchronization: 

N 

ip k = (iu k - y)ip k + 2J Qjkipj, (4) 

where ttjk is the coupling constant of this linear model and 7 is the decay constant to be tuned to 
bring the amplitude of tp k — a complex variable — to a steady state (details given below). As we 
will later show, the argument of ip k corresponds to 9 k in eq. (QQ). For simplicity, we assume 7 > 0, 
and uj k is real. In this reformulation we consider global coupling because it is convenient and will 
permit easy comparison of our results with those of previously studied mean-field models, but this 



approach lends itself equally well to other linear coupling schemes llllh . 
The linear model (eq. H]) has the simple solution 

N 

$=J2 a jVje Xj \ (5) 

3=1 

where aj are constants determined by the initial conditions, and Vj and Xj are the eigenvectors and 
eigenvalues associated with the matrix defined by the RHS of eq. §4§. The synchronizing behavior 
of the system in the long-time limit is dependent on the eigenvalue(s) with the greatest real part. 
Thus in this analysis we will adopt the convention of ordering all eigenvalues by their real part, 
from Ai (least) to A a? (greatest). Distinct eigenvalues with the same real part shall arbitrarily be 
assigned consecutive subscripts within the larger sequence. In our discussion we will assume X N 
is not degenerate. 

We tune 7 so that 3?[Ajv] = 0. If 3?[Ajv-i] 7^ then the solution becomes 

lim if = a N v N e iuJrt (6) 

t^oo 

where the collective frequency of the fully locked state, cu r , is given by 

u r = —iX^ = |A/v| (7) 

and is related to the collective phase of the locked state by <p — w r t. As a result, r will tend to 
a steady-state value between and 1, which (in a finite- N system) indicates full locking, where 
the entire oscillator population is locked to one particular frequency. The ensuing analysis will 
focus on the transition to this fully locked state (the full locking transition). For finite- N systems, 



the transition may be from partial locking — where there are subpopulations locked to different 
frequencies — to full locking, or from incoherence to full locking. 

One can calculate the properties of the steady-state synchronization phase, where only one 
eigenvector remains in the long-time limit. To find an expression for r one must determine the 
eigenvector associated with X N . Hereafter, for simplicity we assume = > unless other- 
wise specified. The general form of the corresponding eigenvector for 1 < j < N, where j is the 
index for the components of v^, is given by 

i(u>N — w r ) — VL/N — 7 



M 3 



(8) 



3 i(uj — u r ) — Q/N — 7 ' 
The general explicit expression for r of a fully phase-locked system with an arbitrary distribution 



of Uk over finite N in the long-time limit is then 
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(9) 



which, it should be stressed, is independent of initial conditions. If we assume the distribution of 
frequencies to be symmetric about uj r , i.e. g(cu r — Uk) = g(uj r + Uk), then we can simplify eq. © 
to arrive at 



1 N 
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(10) 



V^/iV + 7^ 

If r goes to a nonzero steady-state value, there is steady-state synchronization of the system. In a 
finite system, the full locking transition takes place where r goes between having and not having a 
steady-state value. This transition point occurs for a given frequency distribution g(ujf.) when the 
following relationship with the critical value of the coupling constant fi c is satisfied: 



%-i(fi c )]=0. 



(ID 



To see the connection between the linear reformulation, eq. ©, and the Kuramoto model, eq. 
O, we perform the nonlinear transformation ipk(t) = Rk(t)e l9k ^ on eq. © to arrive at 



N R 

$k = uj k + sm(9j - 9 k ) 



n 



N 



R k (t) = -jR k + ^J2 R i cos ( 9 i ~ °^ 



(12) 



(13) 



By tuning 7 so that 3?(Aat) = we can force each Rk to go to a steady state, namely 




(14) 



If all Rk go to a steady state for large times, then eq. (fT2j) becomes eq. (0Q). Therefore, this linear 
version maps onto eq. (Q~|) with the effective coupling constant 



It is important to note that Kj k is independent of initial conditions as cancels out, and that the 
mapping holds only in the regime in which there is steady-state synchronization (where 3ft [A jv-i] 7^ 
0). With this, eq. (fT2j) can be rewritten as 



In other words, by introducing the amplitude R k properly constrained by the decay constant 7, 
we would be able to perform a nonlinear transformation of the Kuramoto model eq. (OQ) (which 
only has an implicit solution in the infinite- iV limit) with an effective coupling constant of Kjk 
into our linear version eq. © that can be solved exactly for any N. One can conceive of much 
more general mappings between eq. © and eq. ©, such as by allowing each oscillator to have 
its own independent 7 — > jk, or by breaking the assumption of uniform coupling in the linear 
reformulation, e.g. f2 — > Vtj k (which does not break the linearity) to accommodate a larger class 
of couplings Kjk. 

IV. EXAMPLES 

To emphasize the mathematical properties of our linear model, we work through three groups 
of examples in detail below. In each, we solve for r using the linear approach and compare the 
result to the Kuramoto solution, eq. ©, in the thermodynamic limit (the regime of validity of the 
Kuramoto solution). 

A. Identical oscillators 

In our first example we consider a system where all characteristic frequencies are equal, i.e. 
uj k = 10. This is a simple system which will exhibit steady-state perfect synchronization regardless 
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(15) 



A' 




(16) 



of initial conditions (as long as 7^ 0). In this case, to solve the linear model (eq. ©), we first 
observe that there are N — 1 degenerate eigenvalues, each equal to icu — 7 — Q/N and one unique 
eigenvalue, 

\ N = ZUJ - 7 + Vl(N - 1)/N, (17) 

that has an associated eigenvector vn = {1, 1, 1, 1}. (Note that the normalization of the 
eigenvector does not matter for the calculation of the order parameter because r is defined such 
that every element has a modulus of 1). In order for each \ip k \ to go to a steady state, we set 
7 = f2(iV — 1)/N (consistent with eq. (TT3T) for a perfectly synchronized state) so that as t — ► 00 
all of the degenerate eigenvectors decay away on a time scale n^j^Q , which for large N reduces 
to ~ 1/fL Assuming Q ^ 0, the solution to the reformulation with oj k = w is then 

\imif =a N v N e iuJt . (18) 

t— >oo 

The synchronization order parameter can now be computed easily, giving the steady- state value 
of r = 1 for long times. Because each element in the sum is normalized by its magnitude, the long- 
time steady-state behavior of r is independent of the initial conditions represented by aj. Since 
Kjk — > 0,/N in the long-time limit in this example, in the infinite-iV limit the linear model's 
solution, eq. (fTOl) . and the Kuramoto solution, eq. (0), give the identical result of r = 1 assuming 
the frequency distribution Uk = oj. It is also instructive to take small perturbations around the 
characteristic frequencies, i.e. uj k = LU + erj k where 0<e< 1 and r] k are symmetrically distributed 
around 0. In this case, to leading order 

^'-4 <i9) 

N on2 



where A 2 is the variance given by A 2 = J2j=i Vj 



B. Bimodal distribution of characteristic frequencies 

For our second example, we look at a system in which both steady-state partial synchronization 
and a full-locking transition occur. Consider a bimodal distribution where u>j = cu — y for 
1 < j < N/2 and Uj = u + f for N/2 < j < N where A > is the width of the distribution 
(we will assume iV is even for simplicity). The solution to this problem maps onto a system of 
two coupled oscillators because of the pecularity of the distribution. 
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Solving the linear reformulation, we see that there are (N — 2)/2 degenerate eigenvalues given 
by i(cu — y ) — — 7 and (N — 2)/2 degenerate eigenvalues given by i(uj + y ) — ^ — 7. The 
final two eigenvalues A at and Aat_i are given by 

Q(N-2) V^ 2 - A 2 
Aat,tv-i = iojq - 7 H — ± (20) 

where the '+' refers to X N , and the '— ' refers to Ajv-i- 

First let us take the situation where Vt < A. Here we shall see that r does not go to a steady- 
state value and therefore that the system will not fully lock and no steady-state synchronization 
occurs. Unlike the previous example, there are not one but two distinct eigenvalues with the largest 
real part, namely 3?[Aat_i] = ^[Xn] = — 7 + ■ So even after setting 7 = n ^ 2) we find 
that 

« WO~ 9 ' l[wo + - o I* 



lim ip = a (A r_ i)W( A r_i)e V / + a N v N e V / . (21) 



It is clear that as long as < A, r will never go to a steady-state value and hence the whole 
system will never fully lock or reach steady-state synchronization. However, independent of all 
initial conditions and assuming f2 > 0, the population of oscillators will split into two perfectly 
synchronized subsets with different collective frequencies given by ^[A^jv-i] =wo± VA 2 -n 2 ^ 

Now let us consider the situation where f2 > A. If we set 7 = ^^T 2 - 1 + ^ n ' 2 ~ A ' 2 (note 
7 — > Cl(N — 1)/N for large fi, consistent with eq. (fT3l) for a perfectly synchronized state) then, 



on a time scale given by 1/ V^ 2 — A 2 , all but one eigenvalue die out over time: 



lim if) = a N v N e iuJot (22) 



where the components of the eigenvector are given by 



= »A-20/iV-7ra_ 
y 13 —iA -20,/N — V^ 2 - A 2 

for 1 < j < N/2 and (v N )j = 1 for N/2 < j < N. Thus r goes to a steady-state value: 



\ 



A\2 



(24) 



2 

for f2 > A. It is important to note that although this equation is true for finite N — unlike the 
Kuramoto solution eq. © — the expression for r is independent of N. 

Since when < A the system does not fully lock and no steady-state synchronization occurs, 
and when f2 > A there is full locking and steady-state synchronization, it is clear that when Q c = 
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A a "first-order" full-locking transition occurs, marked by steady-state partial synchronization 
beyond r c = 1/ v2. Using eq. (|24|) we can deduce a square-root scaling law of the order parameter 
near the full-locking transition (which in this case is also the transition into a steady-state self- 
synchronizing system), i.e. 

i / n-n c 

r - r ^2\j^r- (25) 

As in the previous example, this bimodal example has the property of Kj k — > £l/N in the 
long time limit, so the linear form maps onto the Kuramoto model (eq. (OQ)) with the simple 
transformation — > K. From the above discussion therefore, if one assumes a bimodal frequency 
distribution, i.e. g(uS) = ^i^Vg) _|_ s(^+a/2) ^ solving eq. ©, which is only valid in the 
infinite- N limit, gives our result, eq. (l24j) . 

C. Continuum limit 

In this final example, we work exclusively in the continuum limit and demonstrate that, even 
within this regime where the Kuramoto solution is valid, the linear reformulation still has the ad- 
vantage of a simpler solution for the order parameter and opens up a new class of exactly solvable 
synchronization models. We first determine the critical point of the locking transition and the syn- 
chronization order parameter given by the linear model for a system of oscillators with a general 
symmetric characteristic frequency distribution and a particular effective coupling. Systems with 
such coupling, which we will consider below, go from completely incoherent to fully locked at 
one point in parameter space, i.e. there are no partially locked states; the partial-locking and full- 



locking transitions merge 111511 . Next, we solve for the critical point and the synchronization order 
parameter for a Lorentzian frequency distribution and a uniform frequency distribution. We then 
look at how the result differs from the traditional uniform coupling solution. 

We begin by solving for r, 7, and the locking transition point (also the synchronization critical 
point), each in terms of f2 and a general symmetric frequency distribution. As N — > 00, eq. (flOl) 
becomes 

LO — UJ r 

7 

which holds for Vt > fi c . The anomalous scaling properties of this solution are discussed in 



dujg{uj) 



1 + 



-1/2 



(26) 



II 1711 . The linear reformulation eq. © in the continuum limit has the same form as the equation 



describing the evolution of the fundamental mode in the stability analysis in 111411 of the incoherent 
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state in the Kuramoto model: 

/oo 
g(u')ij)(u' , t)du', 
-oo 

which has an effective coupling in the continuum Kuramoto model of 

K(u,u/) = tt, 



(27) 



(28) 



(u - uj r ) 2 + 7 2 ' 

Following the analysis in the spectrum of the linear operator of the RHS of eq. (1271) comprises 
an eigenvalue \ N at the origin of the complex plane and a continuous line of eigenvalues along 
—7 of the real axis. If Q > Cl c , 7 > and Ajv lies apart from the other eigenvalues, so there is full 
locking and synchronization in the long-time limit. However if fi < $7 C , then 7 = 0, A a? merges 
with the continuum along the imaginary axis, and contribution from every frequency remains in the 
long-time limit, which results in r = as all oscillators are "drifting". This implies a second-order 
phase transition. In this case, the partial-locking transition (often referred to as the synchronization 
phase transition) and the full-locking transition occur at the same point. Furthermore, the coupling 
described by eq. (1281) is such that, at the onset of the transition, there is a uniform distribution of 
phases from to 2n across the oscillator population (i.e. r — > as Q, — > fi c ). We can determine 7 
using the self-consistency equation Ill4ll 



n 



-duj. 



1UJ 



(29) 



At most one solution for 7 exists and 7 > Ill2h . Therefore the system will reach steady-state 



synchronization for 7 > 0, and the critical point Vt c occurs as 7 — > + . This can easily be shown 
to be 

a = • (30) 



7T(/(0)" 

Now, considering a Lorentzian distribution of frequencies about u r , i.e. g(cu — cu r ) = 
7 r[A :2 +(t-^ pj ' f rom e£ l- <l30l) we obtain fi c = A. Integrating eq. (T29l) gives 7; or = — A. So 
where f2 > Q c , 



t'lc 



cos 



1 



1 -1/2 



(31) 



77 \fi — fi c 

in the long-time limit. Similarly for a uniform distribution of frequencies about u r , i.e. g(u—u r ) = 
-^r for I a; — a; r | < 7rA/2 and otherwise, the same steps give us tt c = A, 7^/ = ^ cot(|^), 
and 



Tunif 



( t,t ^ 2"^ ^ Slll!l ' 



7T fi c 

tiU ' 1 2Y1 



(32) 
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in the steady state for > fi c . (It is interesting to note that r uni f > ri or for the entire parameter 
regime.) By contrast, the result one obtains from the Kuramoto solution eq. © for a Lorentzian 
distribution of frequencies (where a partial population of drifting oscillators can occur) is r = 
a/1 — 2A/K Q3j, |4(] for K > K c = 2A and r = otherwise; and for a uniform distribution of 



frequencies one cannot find an explicit solution. 



V. CLOSING COMMENTS 



In closing, two points deserve mention. First, although the dynamics of the linear model are 
in principle completely solvable, they are not the same as the dynamics of the original Kuramoto 
model since the mapping between eq. © and eq. (OQ) only formally holds when all Rj go to a 
steady state. Second, while in this paper we have restricted ourselves only to global coupling, the 
same analysis should be applicable to any linear coupling scheme, including local or asymmetric, 
as long as one is able to determine the largest eigenvalues and eigenvectors of the RHS of eq. ©. 
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